A first principles simulation of rigid water 
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We present the results of Car-Parrinello (CP) simulations of water at ambient conditions and under 
pressure, using a rigid molecule approximation. Throughout our calculations, water molecules were 
maintained at a fixed intramolecular geometry corresponding to the average structure obtained 
in fully unconstrained simulations. This allows us to use larger time steps than those adopted 
in ordinary CP simulations of water, and thus to access longer time scales. In the absence of 
chemical reactions or dissociation effects, these calculations open the way to ah initio simulations of 
aqueous solutions that require timescales substantially longer than presently feasible (e.g. simulations 
of hydrophobic solvation). Our results show that structural properties and difi'usion coefficients 
obtained with a rigid model are in better agreement with experiment than those determined with 
fully ffexible simulations. Possible reasons responsible for this improved agreement are discussed. 



I. INTRODUCTION 

The importance of water in many areas of science has 
motivated a large number of experimental and theoretical 
investigations of the liquid. However, it is only recently 
that x-ray and neutron diffraction measurements have 
come to an overall agreement for properties such as the 
structure of water at ambient conditioned'^. In addition, 
several details of the water structure remain the subject 
of debate, and many dynamical properties of the liquid 
are not yet well understood. 

Theoretical models have played an important role in 
the interpretation of experimental measurements and in 
understanding the physical properties of water'^. Over 
the last thirty years, rather accurate empirical force fields 
have been developed, which can reproduce not only the 
structure but also many dynamical properties of the 
liquid^^ig,. Although empirical models work well for pure 
water under ambient conditions, they are usually diffi- 
cult to generalize to complex solutions or thermodynamic 
states far from ambient conditions. For example, the ma- 
jority of empirical water models that are in current use 
employ potentials that do not change depending on the 
environment. 

In recent years, it has become possible to simulate the 
properties of a liquid entirely from first principles, with- 
out having to resort to fitted potentials. This is due in 
large part to the development of the Car-Parrinello (CP) 
method^ along with the continual increase in high per- 
formance computing resources. Although rather accurate 
and with the potential of being a truly predictive tool, 
CP simulations are much more computationally intensive 
than classical simulations. 

The case of water is particularly demanding for CP 
simulations. The ionic vibrational spectrum of the liquid 
exhibits high frequency modes, i.e. 0-H stretch (3200 to 
3600 cm~d) and H-O-H bending modes (^1600 cm^^). 
Therefore, in order to avoid a coupling between the ionic 



and electronic degrees of freedom, which could cause se- 
vere inaccuracies in a CP simulation, a relatively small 
fictitious electronic mass ( ^ ~ 400 a.u. for protonated 
water) needs to be used*. In turn, the use of small values 
of /i, together with the high kinetic energy cutoff required 
to describe the oxygen pseudopotential in a plane wave 
description, necessitates the use of small integration time 
steps. The time step may need to be as small as 0.08 
fs, which is approximately ten times smaller than what 
is often used in classical MD simulations with empirical 
inter- atomic potentials. This poses a severe restriction on 
the time scales that can be accessed in CP simulations of 
water. 

We note that when using Born-Oppenhcimcr (BO) dy- 
namics (where the total energy of the system is mini- 
mized at each ionic step), one can use larger time steps 
than in CP simulations, since electronic degrees of free- 
dom are not propagated at the same time as ionic coordi- 
nates. However, the accuracy required to reduce system- 
atic errors on the ionic forces so as to have conservative 
dynamics is such that large number of iterations are usu- 
ally necessary to minimize the Kohn-Sham energy at each 
ionic step. Therefore, the gain in efficiency obtained with 
a large time step is more than counter-balanced by the 
computational time requirement for total energy mini- 
mizations. 

It is interesting to note that the problem of integrat- 
ing fast vibrational modes in simulations of liquid water 
has also been encountered in classical MD simulations, 
where the most common approach has been to completely 
eliminate the high-frequency intra-molecular motion by 
using bond length and angle constraints2ii£. Based on 
this approach, a variety of classical water potentials, e.g. 
the TIP series^, are capable of accurately reproducing 
many of the interesting properties of water. In particu- 
lar, results obtained with the rigid water TIP5P potential 
are in very good agreement with a variety of experimen- 
tal measurements such as the structure, the tempera- 
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TABLE I: Details of the simulations. 



Name Geometry p (g/cc) T (K) dt (fs) /i (a.u.) time (ps) 



A 


Rigid 


0.997 


326 


0.07 


765 


16.8 


B 


Rigid 


0.997 


315 


0.24 


1100 


24.5 


C 


Flexible 


0.997 


291 


0.07 


340 


19.8 


D 


Rigid 


1.570 


603 


0.07 


765 


18.6 


E 


Flexible 


1.570 


600 


0.05 


340 


3.0 



ture of maximum density, diffusion, as well as dielectric 
propertiesS*ii*iS. 

With the aim of investigating how to increase the in- 
tegration time step in CP simulations of water and thus 
access larger time scales, we have carried out calcula- 
tions using a rigid water approximation. In this paper, 
we present the results of these simulations and compare 
them to those obtained with flexible water molecules (i.e. 
without imposing any constraints on the geometry of the 
molecules in the liquid), and we discuss the effect of a 
rigid model on the structural properties of water at am- 
bient and high pressure conditions. Our results show 
that an ab initio rigid water model yields faster diffusion 
and radial distribution functions which are less struc- 
tured than those found with a flexible model. Overall, 
the properties computed with the rigid model are in bet- 
ter agreement with experiment than those determined 
with a flexible model. Possible reasons for this improved 
agreement are discussed. In addition, we present a lo- 
calized orbital analysis of the trajectories obtained with 
both a rigid and flexible water model, and we demon- 
strate that the large dipole moment changes in going from 
the gas to the liquid phase are not significantly altered 
by the rigid water approximation. The use of a rigid wa- 
ter model in ah initio simulations opens the way to much 
longer simulations of solutions where chemical reactions 
and dissociation effects do not occur. 



II. METHODS 

In order to examine how the structural and dynami- 
cal properties of water are altered by a rigid water ap- 
proximation, we have performed a series of first princi- 
ple molecular dynamics simulationsi^ of water with and 
without intramolecular bond and angle constraints under 
ambient and high pressure and temperature conditions. 
The simulations consist of 54 water molecules in a peri- 
odically repeated cubic cell with a lateral dimension of 
either 11.74 A or 10.10 A, which correspond to densities 
of 1.00 g/cc and 1.57 g/cc, respectively. At each density, 
we have compared simulations where the intra-molecular 
geometry of the water molecules are rigid to those where 
the geometries of the water molecules are fully flexible. 

In each simulation, the electronic structure was de- 
scribed within density functional theory (DFT)i^4i^ with 
the PBE generalized gradient approximation^^. The va- 
lence wavefunctions and charge density were expanded 



in a plane wave basis, which was truncated in reciprocal 
space at 85 and 340 Ry, respectively. Norm-conserving 
pseudopotentials of the Hamman type were used to de- 
scribe valence-core interactionai2ii&. 

The simulations were performed with the CP 
technique^, which is based on the use of a Lagrangian 
that couples together the system's electronic and ionic 
degrees of freedom, 
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In Eq. n (r) are the Kohn-Sham orbitals, /i is a fic- 
titious mass parameter used to evolve the electronic de- 
grees of freedom in time. Mi are ion masses, Eks is the 
Kohn-Sham energy, fi are occupation numbers, and Ay 
are Lagrange multipliers, used to impose the orthonor- 
mality constraint ^'^*'^j — 5ij. The equations of motion 
derived from Eq. ^are. 
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Central to the CP method is the introduction of the 
fictitious mass parameter, /i, which enables the electronic 
and ionic trajectories to be propagated simultaneously at 
each time step. The particular value of /x is chosen so 
that the dynamics of the electronic degrees of freedom 
occurs on a time scale that is much faster than, and thus 
decoupled from, the ion dynamics. 

For a given system, an appropriate value of ^ can be 
determined by approximating the dynamics of the or- 
bitals generated by Eq.|2]as a superposition of oscillators 
with frequencies 
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where and Cj are Kohn-Sham eigenvalues of occu- 
pied and unoccupied states, respectively. The lowest fre- 
quency obtained from Eq. ^occurs when is the highest 
occupied state (HOMO) and is the lowest unoccupied 
state (LUMO). As discussed in Ref. HI, although Eq. H 
is a rather crude approximation of the true orbital dy- 
namics, it still provides a useful estimate for selecting an 
appropriate value of 

In the case of water, the highest ionic frequency 
is due to 0-H stretching modes, which are at ~3500 
cm~^ within density functional theory. In addition, the 
HOMO-LUMO gap for water within density functional 
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theory is approximately 4.6 eV—. As discussed in Ref. |^ 
this means that values of /z ~340 au are needed to ensure 
a clear separation between the ionic and electronic de- 
grees of freedom in CP simulations of flexible, protonated 
water. However, when the rigid water approximation is 
used, much larger values of fi can be used because the 
bond distance and angle constraints suppress the high 
frequency ionic motions. The various values of ^ that 
we have used for both the rigid and the flexible water 
simulations are listed in Table |I| 

The initial starting configurations of the simulations 
were generated from previous simulations of flexible 
water— by scaling the intramolecular 0-H distances and 
H-O-H angles, while leaving the spatial orientations of 
the molecules unchanged. The particular values for the 
geometry, 0.9926 A for the 0-H distance and 104.6° for 
the H-O-H angle, were chosen to be the same as the peak 
values for the intra-molecular 0-H distance and H-0- 
H angular distributions obtained in the previous flexible 
water simulations^^. These bond distances and angles are 
similar to the ones used in the well-known TIP classical 
water models (0.9572 A, 104.52°)^. In order to maintain 
the distance and angle constraints at every iteration, the 
SHAKE algorithm was used22*2i. 

In simulations A and B (Tabled, the samples were 
initially heated to a temperature of 450 K and cooled to 
300 K for a period of 2 ps. In simulation D, the system 
was heated to a temperature of 900 K and cooled to 600 K 
over 4 ps. The thermostats were then removed and the 
simulations were run under constant energy conditions 
for the times listed in Table 

Maximally localized Wannier functions (MLWF)^^i^^, 
analogous to the orbitals obtained by the Boys localiza- 
tion procedure2i, were used to determine how the rigid 
water approximation may influence the electronic struc- 
ture of the water molecules. By using a recently proposed 
joint approximate diagonalization scheme, we were able 
to compute the MLWFs "on-the-fly" for a large number 
of configurations^^. Following the procedure proposed by 
Silvestreli et al?- , the centers of the MLFWs have been 
used to compute an approximate dipole moment of each 
of the water molecules in the simulation. 



III. RESULTS 

The primary effect of the rigid water approximation 
is best illustrated by comparing the vibrational spectra 
obtained from a rigid and a flexible water simulation. In 
Fig. ^ the power spectra collected from simulations A 
and C are shown. The spectra were computed directly 
from the velocity autocorrelation function with the maxi- 
mum entropy method. The two highest frequency modes 
in the flexible water simulation C are completely absent 
in simulation A. These modes at 1600 cm~^ and 3150 
cm~^ correspond to intra-molecular H-O-H bending and 
0-H stretching, respectively. By suppressing the high 
frequency modes large values of ^ can be used, which 
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FIG. 1: The computed power spectrum of rigid and flexible 
liquid water (see text). The dashed line corresponds to the 
rigid simulation A and the solid line to the flexible simulation 
C. 

in a fully flexible simulation would normally lead to a 
coupling of the electronic and ionic degrees of freedom, 
and in turn to a breakdown in the adiabiticity of the CP 
dynamics® . 

In addition to the obvious absence of high frequency 
modes in the rigid water simulation, there are small dif- 
ferences in the lower frequency range of the spectra be- 
low ~1000 cm~^. For example, the librational modes are 
slightly shifted to lower frequencies by approximately 60 
wavenumbers in the rigid water simulation A. Overall, 
the differences in the low frequency region of the rigid 
and flexible water simulations are small. 

In Fig. 121 the oxygen-oxygen radial distribution func- 
tions (RDFs) obtained from simulations A, B and C are 
shown. Also displayed in Fig.|21 is the latest experimental 
oxygen-oxygen RDF obtained by an analysis of neutron 
diffraction data^, which is nearly identical to the RDF de- 
termined from a recent x-ray diffraction study of wateri. 
The small differences between simulations A and B indi- 
cate that the choice of dt=0.Q7 or 0.24 fs has little effect 
on the oxygen-oxygen RDF when the rigid water approx- 
imation is used. The differences in the height of the first 
peak and the first minimum in the RDFs are most likely 
due to the ~10 K higher average temperature of simu- 
lation A over simulation B (see Table P). However, as 
discussed in Ref. 8, for the simulation times used here 
(15-25 ps), differences of this magnitude are within the 
expected error bars. In Fig. El we also show the oxygen- 
oxygen RDF obtained from simulation C with flexible 
water molecules. The differences between the rigid and 
the flexible water simulations are significant. The first 
peak height is approximately 0.4 higher in simulation C 
than in simulations A and B, and is shifted inward by 
0.03 A. In general, the RDF peak heights and depths are 
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FIG. 2: The oxygen-oxygen radial distribution function of 
water at ambient conditions. The dashed and dotted lines 
correspond to the rigid water simulations A and B, respec- 
tively, the solid line to the flexible water simulation C, and 
the dot-dashed line to experiment!-. 



FIG. 4: The hydrogen-hydrogen radial distribution function 
of water at ambient conditions. The dashed and dotted lines 
correspond to the rigid water simulations A and B, respec- 
tively, the solid line to the flexible water simulation C, and 
the dot-dashed line to experiment-. 
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FIG. 3: The oxygen-hydrogen radial distribution function of 
water at ambient conditions. The dashed and dotted lines 
correspond to the rigid water simulations A and B, respec- 
tively, the solid line to the flexible water simulation C, and 
the dot-dashed line to experiment. 



decreased in the rigid water simulations, showing that 
the rigid water approximation causes an overall decrease 
in the structure of the liquid. The same trend in terms of 
the peak heights and depths has been observed in flexi- 
ble and rigid water simulations by using the SPC classical 
interaction potentiai^S. 

The tendency of the rigid water approximation to de- 
crease the structure of water can also be seen in the 



oxygen-hydrogen and the hydrogen-hydrogen RDFs, as 
shown in Figs. |3| and 01 It is interesting to note that 
in simulation C the second peak in the oxygen-hydrogen 
RDF at r'^2.8 A is larger in height than the third peak 
at r~3.2 A, whereas in simulations A and B, as well as 
in the experimental measurement?;, the second peak is 
smaller than the corresponding third peak. Since the 
second oxygen-hydrogen peak coincides with pairs of hy- 
drogen bonded water molecules in the liquid, the relative 
heights of the second and third peaks have been used as 
an indicator of the amount of hydrogen bonding in the 
liquidSa. 

Overall, Figs. |21 to 0] suggest that the main effect of 
the rigid water approximation is to remove a large frac- 
tion of the over-structure that is characteristic of flex- 
ible water as described within DFT both by CP® and 
BO molecular dynamics''^-, and in general, the rigid wa- 
ter distribution functions appear to be in better agree- 
ment with experiment. There are several possible reasons 
for this improved agreement. For example, it is possi- 
ble that the rigid water approximation mimics some of 
the effects coming from the quantum nature of protons, 
which are explicitly ignored in simulations that allow for 
intramolecular flexibility with classical dynamicE'52i22iS. 
More specifically, in Ref. it was pointed out that at 
a temperature of 300 K, ksT ~ 200 cm~^, whereas the 
high frequency intramolecular modes in water range from 
^1000 to 3500 cni^^. In the quantum system the amount 
of thermal energy available to excite vibrational modes 
is much smaller than the lowest possible intramolecular 
vibrational excitations, 

hiu > fcsT. (5) 
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This indicates that the real quantum system will essen- 
tially be restricted to its vibrational ground state at 300 
K. Therefore, the quantum system may be better de- 
scribed by a classical rigid water model than by a classi- 
cal flexible model, despite the fact that in the quantum 
system the protons become delocalized compared to the 
classical case. 

The effect of the proton quantum motion can be con- 
sidered with path-integral (PI) methods^i. To date, all 
PI simulations with empirical interaction potentials have 
found an overall softening in the RDFs of water at am- 
bient conditionsiii^Si^Si^Si^. As pointed out in Ref. 
these structural changes are quite similar to a 50 K in- 
crease in the simulation temperature. However, it is not 
clear if the lack of quantum effects can fully account for 
all of the overstructure that appears to be characteristic 
in the DFT/GGA treatment of water. In addition, a re- 
cent PI-DFT simulations^ has found the surprising result 
that quantum effects may enhance hydrogen bonding in 
water. 

In order to examine how the dynamical properties of 
water are affected by the rigid water approximation, we 
have estimated the diffusion coefficients, D, for the sim- 
ulations of rigid and flexible water. To compute the 
self diffusion coefficient, the mean square displacement 
(MSD) of the oxygen atoms in the simulations were 
tracked as a function of the simulation time. The sta- 
tistical sampling of the simulation data was improved by 
defining multiple starting configurations separated by 4 
fs. The slope of the resulting MSD in the range of 1 to 10 
ps was then determined and used to estimate D according 
to the Einstein relation24, 

6Di=^lnn |(|r,(t)-r,(0)|^). (6) 

The diffusion coefficient can also be estimated via in- 
tegration of the oxygen velocity autocorrelation function 

as 

^2 = i^°°rfi(vW-v(0)), (7) 

which for the simulation times used here yields nearly 
identical results as Eq. ^ The measured diffusion co- 
efficients for the different simulations as computed by 
Eqs. inland [3 are reported in Table UTI For the rigid water 
simulations A and B, the diffusion coefficients are ap- 
proximately 6.2 and 3.2 times larger, respectively, than 
the flexible water simulation C. Although simulations A 
and B were both performed at higher temperatures than 
simulation C, the observed increases in the diffusion coef- 
ficients are outside the increase (1.4 to 1.7 fold) expected 
for a 15 to 26° increase in temperature over ambient, 
based on experimental data™. Given that the rigid wa- 
ter approximation results in a less structured liquid than 
the flexible one, it is not surprising that it also leads to 
faster diffusion, which is in closer agreement with the 
experimental measurement of 2.4 x 10~^ cm^ / a^^^. 



TABLE II; Summary of computed diffusion coefficients. Di 
corresponds to the diffusion coefficient as estimated by the 
Einstein relation and D2 is from integrating the velocity au- 
tocorrelation function. 



Simulation 


Di (10~^cmVs) 


D2 (IQ-^cmVs) 


A 


2.5 


2.4 


B 


1.3 


1.2 


C 


0.3 


0.4 


D 


2.5 


2.5 




r(A) 



FIG. 5: The oxygen-MLWF center radial distribution func- 
tion for rigid and flexible water. The dotted line corresponds 
to the rigid simulation B and the solid line to the flexible sim- 
ulation C. The dashed vertical lines represent the location of 
MLWFs for an isolated water molecule. 



In addition to possible changes in structure and dy- 
namics, we have examined how the rigid molecule ap- 
proximation changes the electronic properties of water. 
In order to do this we have performed a localized or- 
bital analysis by computing the MLWFs for a series of 
well-separated snapshots from the rigid water simulation 
B and the flexible water simulation C. Within the pseu- 
dopotential approximation, there are four doubly occu- 
pied MLWFs around each of the water molecules in the 
simulations. Two of the MLWFs are localized on the 
oxygen-hydrogen covalent bonds, and the other two are 
localized on the lone-pair locations of the oxygen atoms. 
Given the large amount of data that the MLWFs repre- 
sent, in the following, we only consider the centers of the 
MLWFs rather than the orbitals themselves. In Fig. |31 
the oxygen-MLWF center RDFs for simulations B and C 
are shown. The RDFs consist of two distinct distribu- 
tions centered at r~0.33 A and r'^0.49 A, which corre- 
spond to lone-pair and covalent bond locations, respec- 
tively. For comparison, the dashed vertical lines in Fig. [51 
represent the locations of the MLWF centers around the 
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FIG. 6: The probability distribution of molecular dipole mo- 
ments of water computed from MLWF centers. The dotted 
line corresponds to the rigid simulation B and the solid line 
to the flexible simulation C. 



FIG. 7: The oxygen-oxygen radial distribution function of 
water at p—1.57 g/cc. The dashed line is from the rigid water 
simulation D, and the solid line is from the flexible water 
simulation E. 



oxygen atom of an isolated gas phase water molecule. 
Surprisingly, the rigid water approximation does not sig- 
nificantly alter the large changes in the MLWF centers 
that are expected when going from an isolated water 
molecule to the liquid state. As can be seen in Fig.El the 
lone pair distributions are shifted away from the oxygen 
atoms by ~0.03 A, and the covalent bond distributions 
are shifted toward the oxygen atoms by ~0.04 A. 

As proposed by Silvestrelh et alm^, an approximate 
dipole moment for each water molecule in the liquid can 
be defined by assigning the total charge of each MLWF 
to a point charge located at its corresponding center. Be- 
cause the MLWF on neighboring water molecules in the 
liquid do not significantly overlap, this provides a less 
ambiguous definition of the molecular dipole moments 
than arbitrarily assigning electron density to individual 
water molecules. As pointed out in Ref. |4l|, dipole mo- 
ments computed in this manner from static configura- 
tions may not be representative of the experimentally 
measured dipole moments in the fluid. However, the 
MLWF dipole moments are useful for examining qualita- 
tive differences in the polarization of water as a function 
of different approximations or of solutes present in the 
liquid. 

In Fig. El the probability distributions of molecular 
dipole moments for rigid and flexible water systems calcu- 
lated from the MLWF centers are shown. Fig.^lindicates 
that the rigid water approximation has a rather small ef- 
fect on distribution of dipole moments in the liquid. In 
particular, the average moment in simulation B is shifted 
to 3.08 Debye as compared to 3.20 Debye in simulation C. 
Apparently, an explicit description of high frequency O- 
H stretch and H-O-H bending modes is not necessary to 
reproduce the broad range of moments that are charac- 



teristic of the liquid. It is also interesting to note that the 
latest experimental estimate based on an analysis of the 
x-ray structure factor of water indicates that the dipole 
moment of water in the liquid is 2.9 Debya^^, which is 
in closer agreement with the rigid water model than the 
flexible simulation. 

The decrease in the average dipole moment obtained 
from the rigid water approximation offers another ex- 
planation for the observed softening of the liquid struc- 
ture. In addition to mimicking quantum effects, it is 
possible that the rigid water approximation to some ex- 
tent corrects for the general tendency of simple GGA- 
based functionals to overestimate the polarizability of 
molecule3^2ii^ilSii&. For example, the static isotropic po- 
larizability of an isolated water molecule is 10.74 au with 
the PBE functional as compared to the experimental 
value of 9.64 aup^. It is interesting to note that hybrid 
DFT functionals, which include some amount of Hartree- 
Fock exchange, appear to significantly improve on the 
polarizability of water. In particular, the average po- 
larizability of the water molecule is 9.78 au with hybrid 
PBEO functional'**'. 

In addition to water under ambient conditions, we have 
also examined how the rigid water approximation affects 
the properties of water under extreme temperatures and 
pressures. In particular, we have performed a simula- 
tion of rigid water (simulation D) and of flexible water 
(simulation E) at a density of 1.57 g/cc and an average 
temperature of ~600 K. These high density and temper- 
ature conditions correspond to a regime where molecular 
dissociation is still considered a rare evenliSSiii. How- 
ever, the pressure ('--^lO GPa) is high enough to cause 
a large increase in the nearest neighbor coordination of 
each water from 4.5 at ambient conditions to nearly 13 at 
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FIG. 8: The oxygen-hydrogen radial distribution function of 
water at p=1.57 g/cc. The dashed line is from the rigid water 
simulation D, and the solid line is from the flexible water 
simulation E. 
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FIG. 9: The hydrogen-hydrogen radial distribution function 
of water at p=1.57 g/cc. The dashed line is from the rigid 
water simulation D, and the solid line is from the flexible 
water simulation E. 



lar. In particular, both the large increase in the number 
of nearest neighbors as well as the stiffness of the first 
peak in goo(r) a function of compression are repro- 
duced by the rigid water model^. As higher densities 
and temperatures are considered, intramolecular disso- 
ciation will become an common event in flexible water 
simulations^^ and the rigid water approximation is ex- 
pected to be inappropriate for the description of the liq- 
uid. 



IV. CONCLUSION 

In summary, we have used a series of Car-Parrinello 
molecular dynamics simulations to examine how the rigid 
water approximation affects the computed properties of 
water in the liquid state at ambient conditions. In agree- 
ment with previous observations based on empirical in- 
teraction potentials'^'^, the rigid water approximation is 
found to cause an overall decrease in structure and an in- 
crease in diffusion of the liquid. These changes result in 
properties that are in better agreement with experimen- 
tal measurements than the corresponding first principles 
simulations with flexible water molecules. At higher tem- 
peratures and densities in a regime where intramolecular 
dissociation is still a rare event, the differences between 
simulations where water molecules are either rigid or flex- 
ible become negligible. 

In addition to an improved structural and dynamical 
description of water, the rigid water model enables 
the use of time steps as large as 0.24 fs within the 
Car-Parrinello scheme (i.e. '^3 times larger than in a 
flexible water simulation). A similar conclusion was 
reached in Ref. 48 for first-principle simulations of a 
cytosine molecule in the gas phase. This represents an 
important advantage for first-principle simulations of 
aqueous solutions where chemical reactions do not occur, 
and opens up the possibility of investigating phenomena 
that take place on a long timescale. For example, 
understanding how water orients around a hydrophobic 
solute may require simulations of the order of 100 to 
200 ps. The rigid water approximation presented here 
may prove to be an accurate and efficient approach for 
describing the interaction between a hydrophobic solute 
and water within a first-principles context. 



high pressure^. The oxygen-oxygen, oxygen-hydrogen 
and hydrogen-hydrogen RDFs for simulations D and E 
are compared in Figs. [7| to El Except for the expected in- 
tramolecular differences due to the constraints, the RDFs 
obtained from simulations D and E are remarkably simi- 
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